Submitted to The Astrophysical Journal 



Gravitational Radiation from Hydrodynamic Turbulence in a 
Differentially Rotating Neutron Star 

A. Melatos 

School of Physics, University of Melbourne, Parkville, VIC 3010, Australia 

amelatosOunimelb . edu . au 
and 
C. Peralta 1 

Max- Planck- Institut fur Gravitationsphysik, Albert- Einstein- Institut, Am Muhlenberg 1, 

D-l 44*76 Golm, Germany 

ABSTRACT 

The mean-square current quadrupole moment associated with vorticity fluc- 
tuations in high-Reynolds-number turbulence in a differentially rotating neutron 
star is calculated analytically, the amplitude and decoherence time of 

the resulting, stochastic gravitational wave signal. The calculation resolves the 
subtle question of whether the signal is dominated by the smallest or largest 
turbulent eddies: for the Kolmogorov-like power spectrum observed in super- 
fluid spherical Couette simulations, the wave strain is controlled by the largest 
eddies, and the decoherence time approximately equals the maximum eddy 
turnover time. For a neutron star with spin frequency v 8 and Rossby num- 
ber Ro, at a distance d from Earth, the root-mean-square wave strain reaches 
^-rms ~ 3 x 10 -24 Ro 3 (y s /30 Hz) 3 (<i/1 kpc) -1 . Ordinary rotation-powered pulsars 
(^s ^ 30 Hz, Ro < 1CT 4 ) are too dim to be detected by the current generation 
of long-baseline interferometers. Millisecond pulsars are brighter; for example, 
an object born recently in a Galactic supernova or accreting near the Eddington 
rate can have v s ~ 1 kHz, Ro > 0.2, and hence /irms ~ 10~ 21 . A cross-correlation 
search can detect such a source in principle, because the signal decoheres over 
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the time-scale r c « 1 x Ro _1 (z/ S /30 Hz)" 1 s, which is adequately sampled by 
existing long-baseline interferometers. Hence hydrodynamic turbulence imposes 
a fundamental noise floor on gravitational wave observations of neutron stars, al- 
though its polluting effect may be muted by partial decoherence in the hectohertz 
band, where current continuous-wave searches are concentrated, for the highest 
frequency (and hence most powerful) sources. This outcome is contingent on the 
exact shape of the turbulent power spectrum, which is modified by buoyancy 
and anisotropic global structures, like stratified boundary layers, in a way that 
is understood incompletely even in laboratory situations. 

Subject headings: gravitational waves — hydrodynamics — stars: neutron - 
stars: rotation 



Introduction 



Shortly after the discovery of radio pulsars, spe culation arose tha t the superfluid interior 



of a differentially rotating neutron star is turbulent (lGreensteinlll970l ) . Since then, the theme 



has resurfaced intermittently d uring the quest to understand pulsar rotational irr egularities 



like glitches and timing noise (lAnderson et al.lll978l ; iTsakadze fc Tsakadzdll980l ). Neutron 



star turbulence can be hydrodynamic, takin g the form of a Kolmogorov-like cascade of 



macroscopic eddies at hig h Reynolds numbers (IPeralta et al.ll2005l . l2006bl ; iMelatos fc Peralta 



20071 : IPeralta et al.l 120081 ) . Complicated vorticity patterns of this sort are observed in ter- 



restrial experiments on spherical Couette flow, which undergo transitions to nonaxisymmet- 



vortices (Belvaev et al 


Il979l: Biihlei Il990l: Junk & Eebersl 


20001: Sha & N 


akabavashi 2001; 


Nakabavashi et al. 2002a bl; 


Nakabavashi & Tsuchida 


2005a 


Peralta et al. 


2006a 


2008f). re- 


laminarization (Nakabavashi & Tsuchida 2005b). or Stewart 


son layer disru 
Neutron star ti 


ption ( 


Hollerbach 


2003: Hollerbach et al. 


2004. 


2006: Wei & Hollerbach 


2008). ] 


irbulence can also 



be quantum mechanica l , comprising a se 



( IGlaberson et al.lll974j: |Jou &: Mongiovi 



2002 



Mastrano fc Melatos 



-sustaining tangle of quantiz ed microscopic vortices 
20041 ; Ijou fc Mongiov]|l2006h . excited by bulk two- 



stream instabilities (lAndersson et aU2OO70 , interfacial two- stream instabilities ( Blaauwgeers et al 



20051 ). or meridional circulation ( Peralta et all2005 



2006b; Melatos fe Peralta 



20071 ). In general, macroscopic and microscopic superfluid turbulence appear to trigger each 



other; it is an u nsolved, chicken-or-egg question as to which comes first (jBarenghi et al.ll2001 
Tsubotall20~09[ ). 



Turbulence powered by differential rotation is axisymmetric when averaged over long 
times but nonaxisymmetric instantaneously. Turbulent flows therefore emit stochastic grav- 
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itational waves. In an incompressible fluid, the waves arise mainly from current quadrupole 
(and higher multipole) source terms. In a compressible fluid, the mass multipoles also con- 
tribu te; indeed, they can domina te, e.g. duri ng post-glitch Ekman pumping in a neutron 
star (Ivan Eysden &: Melatod 120081 ) . Recently, iPeralta et al.l (j2006al ) pointed out that there 
exists a fundamental theoretical uncertainty regarding the shape and strength of the grav- 
itational wave signal emitted by hydrodynamic turbulence. The mechanical stress-energy 
in Kolmogorov-like turbulence is contained mostly in large eddies near the stirring scale. 
Naively, therefore, one might expect the gravitational wave signal to look like a 'dirty si- 
nusoid', which reflects circulation on the largest scales and decoheres in approximately one 
rotation period. However, the instantaneous wave strain is proportional to the second time- 
derivative of the stress-energy tensor, and this quantity is greatest for small eddies near the 
dissipation scale, which turn over most quickly. If the latter effect dominates, one might 
expect the signal to resemble white noise. Of course, large eddies match better to low-order 
multipoles than small eddies, and low-order multipo les typically dominate the gravitational 
wave strain far from the source (jWassermanl 120091 ) . A careful calculation is therefore re- 
quired to select between the various possibilities and reliably estimate the detectability of 
the signal. 



In this paper, we undertake such a calculation by combining the formalism of lKosowsky et al 



(120021 ) and lGogoberidze et al.l (120071 ). developed to calculate the gravitational radiation from 
a turbulent, first-o rder phase transition in the early Universe, together with the formalism of 
Wassermanl (120091 ). developed to calculate the gravitational radiation from nonaxisymmetric 
vorticity fluctuations in neutron stars, e.g. due to clusters of quantized superfluid vortices. 
In £j2j we analyze global hydrodynamic simulations of incompressible, shear-driven neutron 
star turbulence to extract the vorticity correlation function, which feeds into the statistics 
of stress-energy fluctuations in the source. We then calculate analytically the current multi- 
pole moments, root-mean-square (RMS) wave strain, and decoherence time of the resulting, 
stochastic gravitational wave signal in §3J The results are applied in §H to estimate the 
detectability of the signal, e.g. with long-baseline interferometers like the Laser Interferom- 
eter Gravitational- Wave Observatory (LIGO), and its polluting effect on continuous- wave 
searches currently under consideration. Hydrodynamic turbulence imposes a fundamental, 
quantifiable noise floor on gravitational wave observations of neutron stars. Astrophysical 
implications, including the rate of gravitational wave braking in different types of neutron 
st briefly canvassed. 
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2. Turbulent vorticity correlations 



Let u(xL,t) be a turbulent vorticity field which fluctuates stochastically with position x 
and retarded time t in the source. The gravitational wave strain generated by the (I, m)-th 
current multipole is proportion al to th e Z-th time derivative of u;(x, t) integrated over the 
source volume, as discussed by iThornd (119801 ) and in §3.21 If the turbulence is stationary, 
the wave strain averages to zero over times longer than the maximum eddy turnover time. 
However, the RMS wave strain is not zero. It is proportional to the autocorrelation function 
(u)i(x.,t)ujj(x',t') s } integrated over the source volume. Here and elsewhere, angular brack- 
ets denote the usual ensemble average. Expressing a>(x, t) in terms of its spatial Fourier 
transform, we can write 



(wi(x,tV*(x',£')) = €ii m e jpq 



d 3 k 



d 3 k' 



,ik'X— ik''X' 



(2vr) 3 J (2tt) 



'^hk' p {v m %t)vl{k\t')) , 



where v(x, t) is the turbulent velocity field satisfying u(pc,t) = curlv(x, t). For stationary 
turbulence, (^(x, t)u*(x.', t')) depends on t and t' only through the combination r = i! — t. 

Global, three-dimensional, numerical simulations of neutron star turbulence driven by 
differential rotation indicate that the turbule nce is approximately isotropic and stationary 
once t he Reynolds numbe r Re exceeds ~ 10 4 (IPeralta et al.ll2005l l2006bl ; iMelatos fc Peralta 
20071 ; IPeralta et al.l 120081 ). Figure [1] presents meridional streamlines of the viscous (left 
panel) and inviscid (right panel) components of a two- comp onent , incom pressible , Hall- 
Vinen-Bekarevich-Khalatnikov (HVBK) superfluid (IHills fc Robertslll977l ) in a differentially 
rotating shell, with dimensionless thickness 5 = 0.3, Rossby number Ro = 0.1, and Reynolds 
number Re = 3 x 10 4 , showing a snapshot of the flow at time t = 4.8Ro~ 1 f2~ 1 , where f2 
denotes th e angular velo c ity of the stellar surface. The simulation parameters are defined 
precisely in IPeralta et al.l (120081 ) . where the numerical algorithm (pseudospectral collocation) 
is also laid out in detail. □ Although the Reynolds number in Figure [T] is ~ 10 7 times 
less than in a realistic neutron star and only just above the threshold for turbulence, it is 
already possible to see that departures from isotropy are limited to the largest scales, i.e. 
~ R*5, where i?* is the stellar radius. This is true even when the shear is stronger; Rossby 
numbers as high as 0.3 have been investigated. Additional pictorial examples appear in 



1 The incompressible a pproximation is acceptable when the turbulent motions are subsonic, as in a neutron 
star JPeralta et al.ll2005h . 



2 Such simulations typically adopt no-slip and no-penetration boundary conditions for the viscous compo- 
nent, perfect slip or no slip for the inviscid component, a Stokes flo w initially, and a mutual friction force of 



the Gorter-Mellink form appropriate for a quant ized vortex tangle ( Gorter fc Mellink 1949t Glaberson et al 
Il974t IPeralta et al. l l2005HAndersson et alJl2007h . 
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Per alt a et al.l (120081 ). Isotropy is expected to increase with Re, as in other turbulent systems, 
but simulations with Re > 10 5 , which would test this claim, are not feasible computationally 
at present. 

The turbulence in Figure [1] is stationary to a good approximation. The streamline 
pattern reorganizes stochastically on the time-scale and the velocity components at a 
fixed point alternate in sign, in such a way that the vorticity averages to the rigid body 
value 2f2 over the long term. This behavior is summarized in Figure [21 The left panel shows 
the meridional streamlines of the viscous HVBK component at four instants in time, each 
separated by 2f2 _1 . The eddies in the flow change noticeably in shape, size, and position 
throughout the simulation (and in the inviscid component, which is not plotted). In the right 
panel, we graph all three vector components of the vorticity at a mid-latitude point versus 
time, after subtracting the rigid body term 2f2. The mean of each component is plotted 
as a horizontal line for comparison. After initial transients die away, i.e. for t > 20f2~ 1 , 
the turbulent vorticity fluctuates stochastically and without bias about its mean value and 
is independent of t. It has zero mean after adjusting for the residual differential vorticity 
2(Afi)z. 

Isotropic, stationary turbulence has a velocity correlation function (v m (k)v „(k')) = 
V(2ir) 3 P mq (k)P(k)5(k — k), where V is the total volume of the system (and drops out at 
the end of the calculation of any physical observable), P mq {k) = 5 mq — k m k q is a projec- 
tion operator perpendicular to the wave vector k, and P(k) is the power spectrum of the 
turbulence, usually a power law. However, for gravitational wave problems, where we are 
interested in temporal fluctuations of the mean-square wave strain (and hence mean-square 
multipole moments), we are obliged to work with unequal time (t ^ t') correlators of the kind 
appearing in (JTJ). As a working hypothesis, in this paper, we assume the standard Kraichna n 



form for high- Re turbulence in three dimensions (jKraichnanl Il959l ; iKosowsky et al.ll2002l ). 
viz. 

(v m (k,t)v* q (k',t')) = V(27r) 3 P mq (k)F(k,t-t')5(k-k') , (2) 

with 

F(k, r) = P(k) exp[-TT7](k) 2 T 2 /A} (3) 

and 

V (k) = (2tt)- 1 / 2 6 1 ^0 3 . (4) 

In (J3J) and (jlj), e is the energy dissipation rate per unit enthalpy (units: m 2 s~ 3 ), and rj(k)^ 1 
is the eddy turnover time at wavenumber k; that is, turbulent motions with wavelength 27i/k 
in any given local region decohere over time intervals longer than ^(k)^ 1 . 

What is P(k) for neutron star turbulence? Alas, there is little one can say with confi- 
dence about this question, given the impossibility of direct measurements and the worrying 
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experience with terrestrial systems, where idiosyncratic features are often imprinted on P(k) 
by boundary layers and other unavoidable global structures even in simple systems (see next 
paragraph). Nevertheless, in order to make progress, we assume that P(k) is a power law, 
P{k) oc k a , and that the power-la w exponent is close to the Kolmogorov va lue for isotropic, 
high- Re turbulence, a = —11/3 (IKraichnanl Il959l ; iGogoberidze et al.l 120071 ). Following the 
standard normalization recipe, we can then write 



P{k) = V 



(5) 



The power law stretches across an inertial range extending from the wavenumber correspond- 
ing to the stirring scale, 

k s = 2n/R„ (6) 
up to the wavenumber corresponding to the viscous dissipation scale, 

k d = [8e/(27z/ 3 )] 1 / 4 , 



(7) 



where v denotes the kinematic viscosity. 



Direct numeric al simulations provide reasonable support for the Kolmogorov scaling 
( IPeralta et al.ll2008l ). In Figure[3l we construct P(k) for the velocity field in Figure[T]from the 



simulation data by summing the squared pseudospectral coe f ficient s |C n / m | 2 corresponding 



to each value of k for the top 10 3 modes (IPeralta et al.ll2005l 120081 ). That is, we plot P(k) 
as a function of fci?*/27T = (n 2 + I 2 + m 2 ) 1 / 2 , where n, I, and m denote radial, latitudinal, 
and toroidal mode indices respectively. The viscous and inviscid HVBK components are 
analyzed in the left and right panels respectively. We first note that the toroidal contribution 
|t^(k, t)| 2 {open circles) dominates P(k), especially at small k, and adheres closely to the 
Kolmogorov scaling (solid line). The poloidal contributions |f r (k, t)\ 2 (open squares) and 
\vg(k., t)\ 2 (open triangles) deviate from the Kolmogorov scaling at small k because isotropy 
breaks down for the largest eddies, as noted before; but, in any case, P(k) is dominated by 
|f</,(k, t)\ 2 at small k. A least-squares fit to P(k) over the inertial range in Figure [3] yields 
a = -3.52±0.35 for the viscous HVBK component (8 < kR*/2n < 38) and a = -3.55±0.25 
for the inviscid HVBK component (16 < kR^jl-n < 53). These results agree surprisingly 
well with the Kolmogorov value a = —11/3, even though the turbulence in Figured] is not 
fully developed, the inertial range stretches over less than one decade in the simulations, and 
there are departures from isotropy at small k. We have verified t hat spectral filtering , whic h 
is implemented in the numerical solver to enhance its stability (IPeralta et al.l 12005k 120081 ) . 
does not warp P(k) significantly for Re = 3 x 10 4 . 

It is important to reiterate that the scaling ([5]) applies to isotropic turbulence in the 
bulk, e.g. grid turbulence far from any walls. It is known from many laboratory experi- 
ments, e.g. in wind tunnels, that P(k) is modified by the presence of anisotropic global 
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structures like boundary laye rs, § to the point where it may not even be a power law 
( ISaddoughi fc VeeravaUil 11994 ) . The modifications are not merely of academic interest; we 
show in §3] that the amplitude of the gravitational wave signal from turbulence is sensitive to 
the form of P(k), e.g. through a. Unfortunately, calculating P(k) accurately is a formidable 
undertaking even in terrestrial situations, where conditions can be controlled, let alone in 
a neutron star. A voluminous literature exists on turbulent cascades in shear and viscous 



boundary layers; see, for example, the review by iRobinsonl (119911 ). Experiments that use 
stereoscopic particle image velocimetry to measure the instantaneous velocity field and hence 
P(k) find that the boundary laye r is populated by co herent structures, like hairpin vortice s 



( IGanapathisubramani et al.ll2003l ) o r wall-wake flows (IPerry fc Marusidll995l;lMarusidl2 001) 



and anomalous Reynolds stresses (jWietrzak fc Lueptowl 11994; IGanapathisubramani et al. 
20031 ). which are inconsist ent with the K olmogorov model (ISaddoughi fc VeeravaUil 1 1994 ). 
The role of stratification (lFernanddll991l). important in a neu tron star, and the interplay 
between shear and buoyant convection (IMoeng fc Sullivanlll994l ) , are also under active inves- 
tigation. Resolving these matters lies far outside the scope of this paper, but it is important 
to recognize them and work towards a better understanding o yer time, e.g. by im proving 
upon the pioneering superfluid spherical Couette simulations of iPeralta et al.l (120081 ). 



Buoyancy supp resses radia l moti on in a neutron star, arguably reducing the turbulence 
to two dimensions. iKraichnanl (119671 ) postulated that two-dimensional turbulence develops 
two inertial ranges: a —5/3 cascade (a = —11/3), which conserves kinetic energy and runs 
to low k, from the stirring scale up to the system scale; and a —3 cascade (a = —5), which 
conserves mean-square vorticity and runs to high k, from the stirring scale down to the dis- 
sipation scale. In a neutron star, the stirring and system scales are approximately the same, 
so the —3 cascade notionally covers a wider k range than the —5/3 cascade. How ever, labo- 



ratory experiments indicate that the situation is more complicated. For example, llida et al. 



( 120091 ) found P(k) oc k~ 4 for horizontal k and P(k) oc k~ 5 for vertical k. Vertical motions 
are suppressed, but the vertical cascade still plays a key role in routing energy t hrough the 
system , by mediating the formation of sheared stacks of pancake-like structures. ISommeria 
(119861 ) showed experimentally that, even when the largest scales are pancake-like, smaller 
scales remain isotropic and Kolmogorov-like, perhaps due to the action of internal gravity 



3 Boundary layers in a neutron star are thin, e.g. Re -1 / 2 (surface Ekman l ayer), Re~ 1,/3 (Stewartso n 
layer tangent to the core), or Re~ 2 / 5 (equatorial Ekman layer) in units of R* ( Peralta &: Melatos 2009). 
Nevertheless, they influence a large volume of fluid by partitioning the flow globally into cells, thereby 
shaping P(k) at low k. Laboratory experiments also reveal transient, ribbon-like streamers, which resemble 
boundary layers, throughout the body of otherwise isotropic Kolmogorov turbulence. The stre amers generate 



anomalous Reynolds stresses an d significant instantaneous departures from isotopy at low k (|Marusidl2001 



Ganapathisubramani et al.ll2003f ). 
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waves (which are weakly damped at the Prandtl numbers Pr ^> 1 expected inside a neu- 
tron star). As a rule, stratified turbulence fossilizes into tw o dimensions when the activity 
parameter I = e/(vN 2 ) drops below « 7 (llida et al.ll2009l ). where N ~ SOOrads -1 is the 
Brun t-Vaisala frequency in a typical neutron star in beta equilibrium (Ivan Eysden fc Melatos 
20081 ). Some of the candidate astrophysical sources considered in §H satisfy this inequality, 



while others do not. We confirm, with the help of order-of- magnitude estimates in §3.3[ that 
the gravitational wave results in §|3] and §H are probably insensitive to the dimensionality of 
the turbulence, at least for the values of a that one might reasonably expect in a neutron 
star. Nonetheless, we emphasize that the question is far from settled and needs to be studied 
more carefully with more sophisticated, compressible numerical simulations. 



3. Gravitational wave signal 
3.1. Current multipole moments 



The gravitational wave strain measured by an observer at a distance d from a source 
can be written in the transverse traceless gauge as a linear combination o f gravitoelectri c 
('mass') and gravitomagnetic ('current') multipoles; see equation (4.3) of iThornd (119800 . 
The latter components take the form 



n jk 



G_ 

Jd 



EE 

1=2 m=-l 



d l S lm (t) 
dt l 



rp&l.lm 

1 jk ' 



(8) 



In equation (JSJ), S (t) denotes the (/, m)-th current multipole moment, written as a function 



of the retarded time t, and Tf£' lm denotes the associated gravitomagnetic tensor spherical 



harmonic, which describes the angular dependence (or beam pattern) of the radiation field. 
For a Newtonian source (slow internal motions, weak internal gravity), like a differentially 
rotating neutron star, the current multipole moments assume the form 



5' 



Irn 



32tt 



(/ + 2)(2/ + i; 



1/2 



(2/ + 1)!! [2(/-l)(/ + l)J 
(i 3 xr w (xx pv) • Y /_1,/m * 



(9) 



32tt 



I + 2 



1/2 



(2/ + 1)!! [21(1 - 1)(Z + 1). 
x J d 3 xr'x • curl(pv)F' m * 



(10) 



where Y' l,lm is a vector spherical harmonic of pure orbital type, Y lm is a scalar spherical 
harmonic, and p(x) is the fluid mass density. Equation (jUJ) corresponds exactly to equa- 
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tion (5.27b) in iThornd (Il980l ). Equation (fTOl) is derived from 
harmonic in terms of gradients of scalar harmonics, viz. 



[1(21 + l)] 1/2 Y l - 



ljm 



rVY lm + l±Y 



i m 



by expressing the vector 



and then integrating by parts (lWassermanll2009l ) . Physically, therefore, the current multipole 
arises from the magnetic component of the velocity field [equation (jUJ)] or, equivalently, from 
the radial component of the vorticity field [equation ffTUl) ]. We only consider incompressible 
turbulence in this paper (see for which the mass multipoles vanish. 

For the sake of simplicity, we take p to be uniform, i . e. p = 3M*/(47ri2^), where M* is the 
total mass of fluid in the stellar interior; cf. IWassermanl (120091 ). A more realistic assumption 
is that p is incompressible (subsonic flow) but stratified gravitationally. However, we avoid 
stratification in this first pass at the problem because we wish to exploit the scalings for 
isotropic Kolmogorov turbulence in §2], which assume uniform p. As noted in §21 the stratified 
problem is much harder. 



3.2. Root-mean-square wave strain 



The mean wave strain at the observer is zero for stationary, isotropic turbulence, as 
discussed in §|2J Therefore, to assess detectability, we compute the autocorrelation function 

C{r) = (h]?(t)h]?(ty) , (12) 

which is positive definite for t = t', reduces to the RMS wave strain /irms = (l^Jkit)] 2 ) 1 ^ 2 
for t = t', and is a function of t and t' only through the combination r = t' — t for stationary 
turbulence. For the sake of simplicity, we compute C{r) for a single multipole (l,m). By 
doing so, we circumvent the following complication: an observer at a particular position, 
doing a realistic detection experiment, sees cross terms in C(t) which mix multipoles together 
(e.g. S 20 S 21 *). The cross terms only vanish when average d over all possib le observer positions 



[by the orthonormality of T? k ' lm ] see equation (2.36) of Thornd (jl980l )]. In this sense, our 



detectability estimates are conservative; the signal from one multipole is obviously a lower 
bound on the total signal. 

The wave strain autocorrelation function is the ensemble average of a product of time 
derivatives of S lm . In the special case of stationary turbulence, where it is always possible 
to (notionally) fix t (t') at some instant during the average and exchange d/dt' (d/dt) with 
d/dr (—d/dr), it is possible to simplify the ensemble average u sing t he following formula 
from turbulence theory, e.g. equation (12) of iGogoberidze et al.l (120071 ): 

'dS lm {t)dS lm {t')\ d 2 



dt 



df 



dr 2 



{S lm (t)S lm (t')) 



(13) 
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Upon combining (jHJ), ffTOl . f|T2|) . and ([TBI , and noting that |T^ 2 ' /m | 2 < 1, we obtain the 
maximum autocorrelation an optimally situated observer can detect: 



C(r) 



•1) 



32tt 



L(2Z+1)!! 



Z + 2 



2Z(Z-l)(Z + l)dr 2 ' 



J2Z pR* 

drr l+3 ' ' • J ^'~' ) 



dr' (rj 



J <e*Jd!>x <y'-(x)x ■ W (x, t)y"»(x')x' ■ „(x', f )') • 



(14) 



Upon substituting ([I])-® into (fill) , to deal with the ensemble average, and performing the 
k' integral over the delta function, we arrive at the expression 



(2^ 



C(r) 



X 



32tt 



(2Z+1)!! 



Z + 2 



2Z(Z -!)(/ + l)dr 2 ' 



j2Z rR* rR* 

1 drr l+3 ' ''■■'< 



cZr' (r')' 



KP(A;)exp[-7r77(A;) 2 r 2 /4] / d 2 x / rf 2 x' 



x / d 2 kY lm *(i)xi exp(ik • x)F' m (xX- exp(-ik • x')(5y - My-) • 



(15) 



Equation ( 1T51) can be simplified analytically in a number of ways. Here, we elect to 
expand the two plane-wave factors in terms of scalar spherical harmonics, viz. 



e ik -i = 4tt^Y 'j -< L (kq) ^ Y LM *(q)Y LM (k) 



(16) 



L=0 



M=-L 



and exploit the orthogonality properties of the spherical harmonics to make progress. In 
( TTBl . ji denotes a spherical Bessel function of the first kind. The six angular integrals (over 
x, x', and k) now factorize easily, and we arrive at the final expression for C(r), 



C(r) 



■1) 



G 2 p 2 

10 r 2 



32tt 



(2Z+1)!! 



Z + 2 



d 



21 



- + 1) dr 21 

dk k 4 VP(k) exp[-7r?7(A;) 2 r 2 /4] 



c lu r 

x I** drr l+3 f ' dr' (r') /+3 1 '"• f ' h ' D "' 1 -.,iS2-2 
o Jo J 

oo oo L L' 

A TslmLM TSlmL'M'* atLML'M' 



X 



EE E E i^wnfti/M^ 



N 



(17) 



L=0 L'=0 M=—L M'=-U 



with 



K, 



ImLM 



LA/*/ 



1-3 



LML'M' 



d z kY LM ik)Y L ' M '*ik)(5 i3 - kikj 



(18) 
(19) 



- 11 - 



One can evaluate Ki and Nij by writing Xi and k\ in terms of Y 10 and y 1 ^ 1 and using 
Clebsch-Gordan coefficients to evaluate the triple products. For example, is nonzero only 
if L — I ± 1 and \M + m| < 1. In what follows, however, we do the integrals directly with 
the help of a symbolic algebra package. 



3.3. Quadrupole 

Let us begin by specializing to the case I = m = 2, where the sums in (JTTJ) are nonzero 
for L = 1,3 and V = 1,3 only. [The same is true for (Z,m) = (2, 1).] Doing the r and r' 
integrals and r derivatives, we find 

C( T ) = %So / dkk- s VP(k)exp[-7i V (k) 2 r 2 /A] 



5625c 10 d 2 . 

x v (k) 4 [12 - 12irr)(k) 2 r 2 + ir 2 r]{k) 4 T 4 }[>ip 51 {kR*) + ^ 53 (kR*)) 2 . (20) 
The integral defined by 

r\x)= [ x dtt a m, (2i) 



o 



when combined with P(k) through the k integral in ( 120]) . expresses mathematically the fact 
that the wave strain is an incoherent sum of randomly phased eddy motions, a subset of which 
are wavenumber-matched to the multipole moment under consideration. It can be shown 
that IC(r)! 1 / 2 is r oughly proportion al to the square root of the mean number of eddies at 



the relevant scale (iWassermanl 120091 ) . 



The maximum RMS wave strain is obtained for r = 0, when there is no turbulent 
dephasing. The k integral in equation (TSUI) runs from the stirring to the dissipation scale, 
i.e. from k s in ([6]) to k^ in ((?]). The ip function is oscillatory, but its envelope grows oc k 3 for 
a = 5. Hence, if r](k) and P{k) are given by (jl]) and ^ respectively, the integrand scales 
oc k~ 3 overall in the limit r — > 0. The RMS wave strain is therefore dominated by motions 
at the stirring scale. Under these circumstances, with A; s k<i, the result is 

0.59G 2 p 2 Ry 

RMS - c i 0rf 2 • 

If the turbulence is powered by differential rotation, the specific (per unit mass) energy input 
per unit time is e = i? 2 (Af2) 3 ( Landau fc Lifshitzlll959l ). □ from which we obtain 



- 5 x ir28 ' life) (tL) 3 (l^r (lo^) 3 ■ < 23 > 



4 This choice of e is conservative. Other possible choices, e.g. e — R^ffiAfl, imply a higher gravitational 
wave strain. 
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The Rossby number is defined as Ro = AQ/Q, as in §2j 

The above results confirm that we can reliably use the quadrupole moment to estimate 
detectability. As /irms is dominated by k B , and the stirring scale is well matched to I = 2, we 
are not missing a lot of power emitted by eddies at large k and showing up at proportionally 
higher /. On the debit side of the ledger, the assumption of isotropy is least valid at the 
stirring scale, as discussed in §|2j Equation (120]) makes it clear why this may ultimately turn 
out to be a serious flaw: /irms depends quite sensitively on the shape of P{k). For example, 
if the exponent of the power spectrum satisfies a > —5/3, the dissipation scale governs /irms, 
not the stirring scale, and equation fl23|) would exhibit different scalings with Q and AQ, 
plus an explicit dependence on the kinematic viscosity of the stellar interior. 

Buoyancy arguably reduces the turbulence to two dimensions by suppressing radial 
motion, as discussed in §21 However, order-of-magnitude estimates suggest that the gravita- 
tional wave signal is insensitive to the dimensionality of the turbulence. Evaluating equations 
(1T71) and (12(1 for a = —5, appropriate for the vorticity-conserving —3 cascade postulated 



by iKraichnanl (119671 ). we find that h RMS changes by less than 20% for k s <C k d . In reality, 



because h RM $ is dominated by motions near k s , where the forward and reverse cas cades cross. 



a effec tively lies closer to —5/3 than to —3, even in the two-dimensional model of iKraichnan 



(119671 ). and the change is even smaller. In all cases, the power-law dependences on e and 
the other quantities in equation (J22"]) are universal, as long as we have a < —5/3; the shape 
of the turbulent spectrum affects just the numerical prefactor in equation (|22p . and even 
then only weakly, unless we have a > —5/3, whereupon the dissipation scale dominates the 
signal. Theory and experiment are united in deeming spectra shallower than P(k) oc &T 11 / 3 
to be extremely rare in nature; certainly, buoyancy acts in the contrary direction to steepen 
P{k). As noted in §21 although neutron stars are strongly stratified, the activity parame- 
ter I = e/(uN 2 ) « 40(Aft /l rad s^fju/l O m 2 s- 1 )- 1 (N/500 rad s' 1 )" 2 falls below the fos- 
silization threshold I w 7 Jlida et al.lExM ) for Aft < 0.6(^/10 m 2 s _1 ) 1/3 (iV/500rad s" 1 ) 2 / 3 



rads -1 . Some of the candidate sources discussed below in §Hdo not satisfy this inequality 
(see Table HI); that is, hydrodynamic turbulence in these sources may be effectively three- 
dimensional despite strong stratification, lending extra support to the results in this section. 



3.4. Higher multipoles 

As a general rule, the current multipole moment depends increasingly strongly on k^ 
as I increases, because higher multipoles match better to small-scale eddies. This effect 
is communicated through the ip function combined with P(k), e.g. as in equation (120"1) . 
Whenever I increases by one, the r derivatives in (|17[) bring down two extra powers of i](k), 
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the r and r' integrals contribute an extra factor of k 2 , and the ip functions contribute new k 
factors too. For I = 3, we find 

o, 7r 3/-r2 ( „2.-.8/3 /-fed 

4ms = 7 2 o3 c i 2rf2 J k dkk-^ir 2 (kR*)+r\kR*)f (24) 
0AlG 2 p 2 R? /3 e s ^ 

, (25) 



c 12 d 2 



with k s <C kd- The integrand in (j24"l) scales oc k~ 5 ^ 3 , so /irms for / = 3 is dominated by 
motions at the stirring scale, just like for I = 2. The ratio of the I = 3 and I = 2 wave strains 
is 0.83(i2*fi/c)Ro, implying that the I = 3 radiation is significantly weaker for any realistic 
neutron star rotating slower than centrifugal breakup. 

For I > 4, one finds that /irms is dominated by motions at the dissipation scale, not 
the stirring scale. Consequently, the dependence of /irms on & and AQ also changes, and 
a new, explicit dependence on the kinematic viscosity appears. For example, for I = 4, one 
obtains /?.rms oc pR^^k l J^d^ 1 oc z/ -1 / 4 , and the ratio of the I = 4 and I = 2 wave strains 
is ~ ( R ii Q/c) 2 Ro' 2 (k r \R f )^ 3 . I n view of the Kolmogorov scaling kd/k s = Re 3//4 [see equation 



(17) of iKosowsky et al.l (120021 )]. we conclude that the I = 4 radiation is weaker than the I = 2 
radiation for Re < (i?*f2/c) _8 Ro -8 , which is always satisfied except possibly in a strongly 
sheared millisecond pulsar, e.g. one born recently in a supernova. 



3.5. Decoherence time 

The decoherence time r c is the time that must elapse before the instantaneous wave 
strains hJjF(t) and hj^(t r ) become statistically uncorrelated. We define it to be the value of 
t = t' — t at which C(r), defined as the ensemble-averaged correlator in (fl2j) . decreases to 
some fixed fraction (say, one quarter) of its maximum (at r = 0). 

The results in §3.31 and §3.41 demonstrate that /irms is generated predominantly by the 
stress-energy in motions at the stirring scale. It is tempting, therefore, to predict that the 
signal decoheres on the ed dy turnover time at th e stirring scale, i.e. the rotation period of 



the star multiplied by Ro (IKosowsky et al.l 120021 ). In fact, the situation is potentially more 



complicated. Looking at equation (12171) . which describes how C(r) decreases with r, we see 
that the integrand contains three terms, which scale oc k~ 3 , k~ 5 ^ 3 , and k~ 1 / 3 . The first, 
positive term is the only one which contributes to the maximum of /irms, i-e. equation (j22|) ; 
its integral is dominated by k s , as discussed above. The second, negative term is the only 
one which can reduce C(r) and cause decoherence. Its integral is also dominated by k s . The 
integral of the third, positive term involves both k s and k& but is dominated by the former 
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(see next paragraph). Hence the final result for r c contains information about both k s and 
k d in general. But, for the I = 2 and I = 3 signals, it turns out that the dependence on k d 
is very weak, unlike for higher I. This result is crucial for the question of detectability, as 
explained in §HJ 

By graphing C (r) numerically in Figure HI we find that it falls to one quarter of its 
maximum over a time comparable to the eddy turnover time at the stirring scale, ^(fc s ) _1 . 
This graphical task is tricky, because the square of the ip ab functions oscillates rapidly for 
k d ^ k s . However, we can obtain an excellent analytic approximation by averaging over 
many cycles of the fast oscillation. Upon writing [ip bl (x) + tp 53 (x)] 2 ~ 25x 6 /2 plus fast 
oscillations oc cos(2x) in the regime x ^> 1, we arrive at the formula 



C{r) 
h 2 

"-RMS 



)V 



cxp 



7ir](k s ) 2 T 2 



+27r 2 7](k s yr i { Erf 



k s 



Erf 



(26) 



where Erf(x) symbolizes the error function. Equation ( |26i) is accurate to ~ 12% across the 
full range of r for the parameter range under consideration. The dependence on k d /k s 3> 1 
in the third term is weak, although this stops being true for the rare case of non-Kolmogorov 
spectra with a > —5/3, where k d dominates (see §3.31) . From (j4]), (jSJ), ((7j), and (|26|) . with 
e = Rl(AQ) 3 , we arrive at the following expression for the decoherence time corresponding 
to the half-strain point: [f| 



0.35r/(fc s )" 1 



26 ms 



10rads~ 



-i 



(27) 
(28) 



The s hear viscosity in a neutron star is given by the neutron-neutron scattering formula 
derived bv lCutler fc Lindbloml Jl987h . viz. v = 10 (p/6 x 10 17 kg m - 3 ) 5 / 4 (T/10 8 K)~ 2 m 2 s -1 , 
where T is the temperature of the stellar interior. From equation ([7]), we get k d 3> k s for 
fiducial values of p and T. In a protoneutro n star, discussed in §4^ the bulk viscosity (oc 
T 8 ) exceeds the shear viscosity (1Sawyerlll989l ). Under these conditions, small compressions 
enhance the turbulent dissipation rate, effectively reducing k d . Numerical simulations of a 
compressible HVBK superfluid (outside the scope of this paper) are needed to quantify this 
effect properly. However, the foregoing gravitational wave calculations continue to hold to 



5 It should be noted that rj(k) does not always follow © near the stirring scale (jGogoberidze et al. 
as discussed in We aim to refine (f!?T)) in future work. 



20071) . 
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a good approximation, provided that the stirring and dissipation scales remain moderately 
well separated [e.g. k d /k s > 5 in equation (|26j) ]. 



4. Detectability and astrophysical implications 

In this paper, we calculate analytically the current multipole moments generated by 
vorticity fluctuations in high-Re turbulence in a differentially rotating neutron star. We 
derive an analytic expression [equation (120]) ] for the wave strain autocorrelation function 
C{t) of the resulting, stochastic gravitational wave signal to leading (quadrupole) order, in 
terms of the turbulent power spectrum and eddy turnover time spectrum, and show that 
C(r) is governed by motions at the stirring scale. From the same equation, we also compute 
the root-mean-square wave strain /irms and the decoherence time r c of the signal and show 
that r c approximately equals the eddy turnover time at the stirring scale (usually much 
longer than the rotation period). Convenient formulas for h RMS and r c , scaled in terms of 
astrophysical parameters, are presented in (|2"3"]) and (]2"g]) respectively. 



We can use these formulas to assess qualitatively whether the stochastic signal can be 
detected by existing and planned long-baseline gravitational wave interferometers. Clearly, 
neutron star turbulence imposes a fundamenta l, un a voida ble, astrophysical noise floor on 



continuous- wave searches for neutron stars, e.g. labbl (j????l ). But how seriously does it pol- 
lute such searches? To answer this question, we note two things. First, most neutron stars 
are either rotating too slowly or with too little shear to cause trouble, according to (|23|) . 
at least at the sensitivities anticipated for Enhanced and Advanced LIGO. Second, for the 
small subset of neutron stars that are potentially powerful stochastic emitters, the signal de- 
coherence time is relatively short; indeed, for AQ ~ Q and the fastest rotators, r c approaches 
(without ever quite reaching) the sampling time of LIGO-like interferometers . Equations 
fl23|) and (128]) imply that r c decreases as /irms increases. 

Tabled] illustrates these conclusions by listing h RMS and r c for several realistic categories 
of neutron star sources. 



1. Protoneutron stars. The stellar angular velocity profile is a key o utput of radia- 



tion (maKneto)hydro dynamics simulations of core-collapse supernovae (jOtt et al.ll2006 



Burrows et al.ll2007t ). All progenitor models tested so far lead to strong diffe rential ro 



tation. Looking for exa mple at the fast e st rot ators in Figures 8, 9, and 15 of lOtt et al. 



(120061 ) or Figure 14 of iBurrows et al.l (120071 ) . captured 0.2 s after core bounce in a 
two-dimensional, unmagnetized explosion, we see that Q decreases gradually from 
k 3 x 10 3 rads _1 at r = 3 km to ~ 1 x 10 3 rads _1 at r = 30 km (enclosed mass 
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w 1.2M ), before dropping steeply to w 20rads~ 1 at r = 300km (enclosed mass 
« 1.7M ). Conservatively, for a hypothetical supernova in the Milky Way, the 
above figures imply /irms = 9 x 10~ 21 and r c = 0.1ms (first line of Tabled]). LIGO 
is capable of detecting such stochastic emission, because r c is greater than the sam- 
pling time of the interferometer (« 60 /is), but the signal partially decoheres. The 
prospects are brighter if the protoneutron star rotates slower and with less shear, and 
the gravitational radiation emanates mainly from the extended envelope (second line 
of Tabled]), so that lower /zrms (5 X 10 -22 ) is traded for higher r c (3ms). Detection 
by Advanced LIGO is t hen possible in principle, e.g. via a cross-correlation search 
(IDhurandhar et al.1 120081 ). The signal persists for > 10 2 s ^> r c before the differen- 
tial rotation dissipates and/or t he protoneutron star spins down , e.g. due to r-modes, 
fallback, or a magnetized wind (Lai et al.ll200ll ; lOtt et al.ll2006l ). 



2. ditchers. Discontinuous spin- up events ('glitches') observed in some rotation-powered 
pulsars are adduced as evidence that neutron stars rotate differentially; the nuclear 
l attice spins down elec tromagnetically, lagging the superfluid due to vortex pinning 



(lAnderson fc Itohlll975l ). The fractional jump in angular velocity ranges from 10 to 



10 4 across the pulsar population. The surprising absence of a 'reservoir effect' (i.e. 
glitch size cx time elapsed since the preceding glitch) in most obj ects implies that the 



observed spin ups are a smal l percentage of the underlying shear (jMelatos et al.l 12008 



Warszawski fc Melatosll2008l ). It is therefore safe to expect Ro > 10 4 in some pulsars, 



although the proportion is hard to quantify. Two examples are given in Table [TJ an 
adolescent, Vela-like pulsar (age ~ 10 4 yr), which spins relatively slowly (~ 10 Hz) 
but undergoes relatively large glitches (fractional jump ~ 10~ 6 ); and a young, Crab- 
like pulsar (age ~ 10 3 yr), which spins faster (~ 30 Hz) but undergoes smaller gliches 
(fractional jump ~ 10~ 8 ). Assuming that the typical glitch resets ~ 0.01% of the 
underlying shear, we find that the decoherence time (0.3 s to 10 s) matches well to 
the LIGO pass band, but the signal is too weak (/irms ^ 10~ 30 ) to be detected by 
interferometers u nder development, or to pollu te continuous- wave searches targeted at 
glitching pulsars ( van Eysden fc Melatos 20081 ) . 



3. Accretors. The hydromagnetic accretion torque acting on a compact star in a mass- 
transfer binary is often comparable to, or greater than, the electromagnetic spin-down 



Strongly magnetized models rotate ~ 10 times slower (IHeger et al 



millisecond protomagnetar engine for long-duration gamma-ray bursts (jBucciantini et al.ll2007t h 



2005 



Burrows et all 120071): cf. 



7 The persistent, stochastic, gravitational wave emission from shear-driven hydrodynamic turbulence is 
unrelated to the putative burst emission from the glitches themselves. 
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torque acting on an isolated pulsar (jHartman et al.ll2008l ) and fluctuates by several per 
cent daily, causing X-ray variability. Accreting objects are therefore likely to rotate 
differentially, with Rossby numbers comparable to, or greater than, those of isolated 
glitching pulsars. Two examples are given in Table [TJ a standard, accreting millisecond 
pulsar (fifth line), and an accreting white dwarf w ith an ONeMg core (sixth line) in a 



system on the verge of accr etion-induced collapse (IBlackman fc Yilll998l ; iDessart et al. 



20071 ; iMetzger et al.l 120081 ). The angular velo city profile of t he wh ite dwarf before 



(and after) collapse is plotted in Figure 12 of IDessart et al.l (120061 ). In both cases, 
^rms ^ 10 -26 falls below the threshold for detection by Advanced LIGO but may 
approach the sensitivity of the next generation of interferometers. 

4. Nearby neutron stars. The neares t radio pulsar discovered to date is PSR J0108— 1431, 
with d = 85 pc (iTauris et al.lll9941). T he nearest millisecond pulsar is PSR J0437— 4715, 
with d = 157 pc (jVerbiest et al.ll2008l ). However, many radio-quiet neutron stars with 
d < 85 pc should reside in the Solar neighborhood; the evidence is both observational 
[e. g. radio-quiet X-ray point sources like the 'Magnificent Seven', some with parallaxes; 



sec 



Popov et al.l (120031 ) and references therein] and theoretical [e. g. population synt hesis 
models predict ~ 10 9 compact objects in the Milky Way; see Kiel et al.l (120081 ) and 
references therein]. If the nearest objects have d ~ 10 pc and rotate reasonably fast, 
they represent promising LIGO candidates. § For example, Table [1] quotes /irms f° r 
two nearby isolated pulsars with Ro = 10 _1 . The millisecond pulsar in particular is a 
bright source, although it suffers from the usual drawback of fast rotators: r c is short. 
Although it is sometimes assumed that millisecond pulsars do not rotate differentially, 
because they are not seen to glitch, it is equally possible that they glitch strongly 



but i nfrequently, because the electromagnetic spin-down torque is weak (IMelatos et al. 



2008|). 



Even before detecting gravitational waves, we can use existing radio timing observations 
of neutron stars to test and constrain the model in this paper. The radiation emitted by 
hydrodynamic turbulence exerts a reaction torque, which fluctuates in sign instantaneously 
but drains rotational kinetic ener gy from the fluid over time. Th e angular velocity of the 
star decreases at the average rate (lThornelll980l ; I Wassermanl 120091 ) 



\n 



GW 



5G 




Ql+lglm 



dt l +l 



(29) 



8 Our inability to measure an ephemeris from radio timing observations does not affect cross-correlation 
searches for the stochastic radiation from neutron star turbulence, although it is a major obstacle to coherent 
continuous- wave searches. 
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10 km 



10 3 rads -1 



AQ 



10 rad s 



-i 



(30) 



rads" 2 .(31) 



Equation fl3~Tl) follows from (1291 via the sequence of steps in §3.1I - §3.3I In the last column of 
Table [U we compute |^gw| for the source categories discussed above. In several instances, 
the fiducial value of |^gw| is close to, or even greater than, the value of \Ct\ measured in 
radio timing experiments. Already, this places constraints on the source parameters, as the 
inequality |^gw| < 1^1 must always hold for any specific object. For example, the oldest 
radio m illisecond pulsars are m e asured to have \Q\ ~ 10~ 15 rads~ 2 [see, for example, Table 
10.1 in iLyne fc Graham-Smithl (119981 )]. Accreting millisecond pulsars in low- mass X-ray 
binaries have \ Cl\ < 10~ 13 rads~ 2 du ring X-ray outbursts and ~ 10~ 16 rads~ 2 between X- 
ray outbursts (lHartman et al.ll2008l ). If these data are representative of the entire millisecond 
pulsar population, then we can start to rule out the existence of large shears like those quoted 
in the fifth and seventh entries of Table [H Indeed, equation (j3~Tl) and the data combine to 
yield a direct, observational upper limit on the rotational shear in any neutron star, viz. 



Ro < 4 x 10~ 2 



10 3 rad s" 1 



-7/8 



\n\ 



lO" 15 rads" 2 



(32) 



Equation (1321) is an important result. It allows the theory in this paper to be falsified, if 
a glitching pulsar is discovered whose fractional angular velocity jumps exceed the right-hand 
side of (1321) . We will investigate this application more thoroughly in a forthcoming article. 
Figure [5] gives a taste of what is possible. In the left panel of Figure El we plot Ro max [i.e. the 
right-hand side of fl32|) ] versus spin-down age f2/(2|$7|) for all ob jects with Rom^ < 1 in the 



rignt-nand side ot lljzp j versus spm-down age \i/{z\\i\) tor all ob jects witn i\o ma * <-. l m tne 
Australia Telescope National Facility (ATNF) Pulsar Catalog jManchester et al.l l2005h . @ 



Clearly, millisecond pulsars with ages > 10 8 yr place strict limits on Ro, with Ro max ~ 10 -2 
in some cases. Furthermore, the pulsars marked with diamonds have been observed to glitch 
at least once. If the theory in this paper is correct, then the largest angular velocity jump 
observed, (AO/Q) g)1Jiax , cannot exceed the underlying shear, AQ/Q, in each of the glitching 
pulsars; in fact, it is expected to be much smaller, given the absence of a reservoir effect. 
As a test, we plot (Af2/f2) g;max /Ro max versus spin-down age in the right panel of Figure 
El In all cases, the plotted ratio is much smaller than unity, consistent with theory. 
Future observational campaigns to simultaneously monitor large groups of pulsars, e.g. with 



9 Objects with Ro max > 1 do not place a useful limit on the shear. 

10 The limits derived from Figure [5] are conservative. In the standard vortex unpinning paradigm, the 
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multibeam radio telescopes like the Square Kilometer Array, will challenge the theory more 
keenly. 

In summary, hydrodynamic turbulence imposes a fundamental noise floor on gravita- 
tional wave observations of neutron stars. Its polluting effect is muted by partial decoherence 
in the hectohertz band, where current continuous-wave searches are concentrated, but only 
for the fastest rotators with the strongest shear. In addition, the mechanism sets a fun- 
damental lower limit on the spin-down rate \Cl\ and hence an observational upper limit on 
the Rossby number Ro, when combined with pulsar timing data. These conclusions hold 
subject to one major caveat, which is discussed at length in ^J2j to wit, that /irms and r c 
depend somewhat on the exact shape of the turbulent power spectrum, P(k). It is known 
from laboratory experiments that P(k) is modified away from its isotropic, Kolmogorov form 
by the presence of anisotropic global structures like turbulent boundary layers, which are 
profoundly difficult to model in terrestrial contexts, let alone in a neutron star. Likewise, 
buoyancy modifies P(k) by creating unequal turbulent cascades parallel and perpendicular 
to the stratification direction, in a fashion that is still not understood completely even in 
controlled laboratory experiments. The order-of-magnitude estimates in §3.31 provide some 
comfort that the gravitational wave results are insensitive to the above considerations, but 
we emphasize that much work (e.g. compressible HVBK simulations) still needs to be done 
to clarify the issue. 

More complete detectability estimates referring to specific search pipelines lie outside 
the scope of this paper. Likewise, we defer calculating the detailed frequency spectrum of 
the stochastic signal, a substantial task. 
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fractional angular velocity jump observed in a large glitch can be related to the internal shear according to 
~ (7 s /J c )(Ar/i?)(AO/n), where I s /I c ~ 10 2 is the ratio of the superfluid and crustal moments of inertia , 
and Ar/R ~ 10~ 6 is the normalized radial distance moved by the unpinned vortices ([Alpar et al.l ll986). 
Arguably, therefore, the theory is challenged if (Af2/Q) g)max /Ro max exceeds ~ 10~ 4 rather than unity, a 
stronger test. 
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Type 


d (kpc) 


n (rads- 1 ) 


Ro 


h-RMS 


r c (S) 


Ogw (rads 2 ) 


Protoneutron star 
















(a) core (30 km) 


10 


2 x 10 3 


1 


9 x 10- 


-21 


1 x 10~ 4 


3 x 10° 


(b) envelope (300 km) 


10 


1 x 10 3 


0.1 


5 x 10" 


-22 


3 x IO" 3 


1 x 10~ 6 


ditching pulsar 
















(a) Vela-like 


1 


1 x 10 2 


io- 2 


5 x 10- 


-31 


3 x IO" 1 


3 x IO" 27 


(b) Crab-like 


1 


2 x 10 2 


10" 4 


4 x 10- 


-36 


1 x 10 1 


4 x 10" 41 


Accretor 
















(a) millisecond pulsar 


1 


3 x 10 3 


io- 2 


1 x 10- 


-26 


9 x IO" 3 


7 x IO" 17 


(b) white dwarf 


0.1 


0.2 


0.1 


2 x 10- 


-28 


1 x 10 1 


3 x IO" 29 


Nearby pulsar 
















(a) fast rotator 


0.01 


3 x 10 3 


0.1 


1 x 10" 


-21 


9 x IO" 4 


7 x IO" 9 


(b) slow rotator 


0.01 


2 x 10 2 


0.1 


4 x 10- 


-25 


1 x 10~ 2 


4 x 10~ 17 



Table 1: Candidate source categories. 



Fig. 1. — Hydro dynamic turbulence in a two-component, incompressible, HVBK superfluid 
in a differentially rotating shell, with dimensionless thickness 5 = 0.3, Rossby number Ro 
Afl/Q = 0.1, and Reynolds number Re = 3 x 10 4 . Parameters are defined in iPeralta et al. 



( 120081 ) . The figure shows a snapshot of the meridional st reamlines of the two components at 
time t = 48fi _1 , taken from the numerical simulations in IPeralta et ali ( 2008 ). The rotation 
axis points vertically. Left panel. Viscous component. Right panel. Inviscid component. 
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Fig. 2. — Stationarity of the turbulence in Figure HJ for the same simulation parameters. 
Left panel. Snapshots of the meridional streamlines of the viscous HVBK component at (a) 
t = 50fi -1 , (b) t = 52fi _1 , (c) t = and (d) t = 56fi _1 . The streamlines reorganize 

stochastically on the time-scale Q' 1 . Right panel. Vorticity at a fixed point (r = 0.85-R*, 
9 = 7r/4, (p = 0) as a function of time (150 < Qt < 250). The spherical polar vorticity 
components are plotted after subtracting the rigid body rotation, viz. u r — 2Q r (solid curve), 
uj0 — 2Qg (dashed curve), and uj^ (dotted curve). The horizontal lines represent the mean 
value of each component over the plotted range. The vector sum of the means equals 2(Afi)z, 
the residual differential vorticity. 
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Fig. 3. — Turbulent power spectrum P{k) (in arbitrary units) for the simulation parameters 
in Figured], viz. 5 = 0.3, Ro = 0.1, and Re = 3 x 10 4 . The power spectrum is proportional to 
Ylm i m \ Cnim\ 2 i where {n, I, m) are mode indices in the pseudospectral expansion; for a given 
wavenumber bin [k, k + dk], the sum runs over all indices satisfying k < {n 2 + I 2 + m 2 ) 1 / 2 < 
k + dk (among the top 10 3 modes). A snapshot of the spectrum at t — 50fi _1 is plotted 
as a function of normalized wavenumber kR*/2ir, with logarithmic binning. It contains 
three contributions, \v r (k, t)\ 2 {open squares), |fe(k, t)\ 2 {open triangles), and |f<^(k, t)\ 2 {open 
circles), whose sum is proportional to P{k). The Kolmogorov scaling {solid line), with 
arbitrary normalization, is overplotted for comparison. Slightly different bins are used for 
the three terms. The steady-state differential rotation, contained in Coio, lies off the scale. 
Left panel. Viscous HVBK component. Right panel. Inviscid HVBK component. 
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Fig. 4. — Wave strain autocorrelation function C(r), defined by ( 1261) . normalized by the 
mean-square wave strain ^r MS , as a function of the time lag r, normalized by the maximum 
eddy turnover time ^{k^ 1 at the stirring scale, for k<±/k s = 10 3 . 
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Fig. 5.— Upper limit on the angular shear in rotation-powered pulsars, assuming gravi- 
tational wave spin down due to hydrodynamic turbulence. Left panel. Maximum Rossby 
number, Ro max , given by the right-hand side of equation (|32|) . versus characteristic age, 
fi/(2|fi|) (in yr), for objects with Ro max > 1 in the ATNF Pulsar Catalog. The points 
marked with diamonds indicate pulsars with a history of glitch activity. Right panel. Max- 
imum observed glitch size, (AO/f2) gjmax , expressed as a fraction of Ro max , plotted versus 
characteristic age, for the points marked with diamonds in the left panel. The ratio should 
not exceed unity for any pulsar. 



